#R SCRIPT TO REPLICATE ANALYSIS AND RESULTS IN 'MEDIA FREEDOM AND THE INSTITUTIONAL UNDERPINNINGS OF POLITICAL KNOWLEDGE', BY MARTIJN SCHOONVELDE
#FOR QUESTIONS, PLEASE CONTACT ME AT MSCHOONVELDE@GMAIL.COM.

rm(list=ls(all=TRUE))

#load required functions in R
library(foreign)
library(arm)
library(ggplot2)
library(coefplot)
library(Zelig)
library(xtable)
library(Hmisc)
library(plyr)
library(stargazer)
library(reshape)

#corstarsl() can be obtained from this url: http://myowelt.blogspot.co.uk/2008/04/beautiful-correlation-tables-in-r.html; make sure to use correct source path when loading into R
source("~/corstarsl.R")

#load replication data; make sure to use correct source path when loading into R
load("~/Replication.RData")

load("/Users/hjms/Documents/Publications/PSRM/ReplicationData/Replication.RData")

#Measurement of the variables is described in their labels; these can be accessed through the describe() function in R
describe(data)

#create vector of variables used in analysis
vars <- c("sophistication", "income","education", "union", "age", "parliamentary", "compvoting1", "partylist1","majoritarian","enp", "freedom01", "herfindahln", "herfindahlt","newspaper","public", "gdp")

#Replicating Table 2: Summary Statistics
stargazer(data[,vars])

#Replicating Table 3: Correlations
xtable(corstarsl(data[,vars]))

#Running Models 1 through 5 in the paper; for description of each model, see paper

model.vars <- c("A1004","sophistication","income01_wi_01","education01_wi_01", "union", "age", "parliamentary", "compvoting", "partylist","majoritarian","enp", "freedom01", "herfindahln", "herfindahlt","newspaper","public", "gdp", "edxfree") 

model.data <-  data[, model.vars]

Model1 <- lmer(sophistication ~ income01_wi_01 + education01_wi_01 + union + age + freedom01 + parliamentary + compvoting + partylist + majoritarian + enp + herfindahln + herfindahlt + newspaper + public + (1| A1004), model.data) 

Model2 <- lmer(sophistication ~ income01_wi_01 + education01_wi_01 + union + age + freedom01 + parliamentary + compvoting + partylist + majoritarian + enp + herfindahln + herfindahlt + newspaper + public + gdp + (1| A1004), model.data)

Model3 <- lmer(sophistication ~ income01_wi_01 + education01_wi_01 + freedom01 + parliamentary + compvoting + partylist + majoritarian + (1| A1004), model.data)

Model4 <- lmer(sophistication ~ income01_wi_01 + education01_wi_01 + freedom01 + parliamentary + compvoting + partylist + majoritarian + gdp + (1| A1004), model.data)

Model5 <- lmer(sophistication ~ income01_wi_01 + education01_wi_01 + freedom01 + parliamentary + compvoting + partylist + majoritarian + gdp + edxfree + (1 + education01_wi_01| A1004), model.data)

#Replicating Figure 1: Bivariate relationship between voter knowledge and media freedom
mean_data <- ddply(data, .(A1004), summarise, freedom01 = mean(freedom01, na.rm = TRUE), sophistication = mean(sophistication, na.rm = TRUE))

p <- ggplot(mean_data, aes(x=freedom01, y=sophistication, label=A1004))
p <- p + geom_point(size = 1) + stat_smooth(method="lm", color="black") + geom_text(hjust=0,vjust=0,size = 2.25) + scale_x_continuous("Media Freedom") + scale_y_continuous("Average Voter Knowledge") + theme_bw()
p <- p + theme(axis.title.x = element_text(face="plain", color="black", size=15, hjust = .5,vjust = 0.1))
p <- p + theme(axis.title.y = element_text(face="plain", color="black", size=15, hjust = .5, vjust = 0.1))
print(p)

#Replicating Figure 2: Conditional Effects
#NB: The code below generates a figure that is slightly different from Figure 2 in the paper but which leads to identical substantive conclusions. Figure 2 contains two small errors in that (i) the standard errors are too large due to a coding mistake, and (ii) the conditional effects are off.

ceplotrvars <- c("A1004", "sophistication","education01_wi_01","freedom01","edxfree","zeta2")

ceplot <- data[,ceplotrvars]
ceplot <- na.omit(ceplot)

attach(ceplot)
ceplot$Betaj <- fixef(Model5)[[3]] + fixef(Model5)[[10]]*freedom01 + zeta2
#The Betaj vector represents the conditional effects (see Equation 2 in paper) with zeta2 denoting the random effect of education within each country.

y <- ddply(ceplot, .(A1004), summarise, freedom = mean(freedom01, na.rm=TRUE), education = mean(education01_wi_01, na.rm=TRUE) , soph = mean(sophistication,na.rm=TRUE), conditional = mean(Betaj), zeta2 = mean(zeta2))

mixedEffect <- melt(as.matrix(ranef(Model5)[[1]]))
mixedEffect$SE <- melt(as.matrix(se.coef(Model5)[[2]]))[, 3]

mixedEffect1 <- mixedEffect
mixedEffect1 <- mixedEffect[32:62,]
levels(mixedEffect1$X2) <- c("Random Intercept", "Random Coefficient")
mixedEffect1 <- cbind(mixedEffect1, y)

mixedEffect1 <- transform(mixedEffect1, X1=reorder(X1, freedom) )

png(file="Figure2.png", width=2212,height=1350,res = 300)

zp3 <- ggplot(mixedEffect1,  
              aes(x = X1, y = conditional,
                  ymin = conditional - 2*SE,
                  ymax = conditional + 2*SE))
zp3 <- zp3 + facet_grid(. ~ X2,panel.margin = unit(1, "lines"))
zp3 <- zp3 + geom_pointrange()
zp3 <- zp3 + coord_flip()
zp3 <- zp3 + theme_bw()
zp3 <- zp3 + ylab(NULL) + xlab(NULL)
zp3 <- zp3 + theme(plot.title = NULL)
zp3 <- zp3 + stat_smooth(aes(x=mixedEffect1$freedom*31, y=mixedEffect1$conditional), method="lm",size=.25, alpha = 0.20,colour=I("#821122"))
zp3 <- zp3 + scale_y_continuous(breaks=seq(-1,10,1/20))
print(zp3)

dev.off()

#regression of conditional effects on media freedom
reg.line <- lm(conditional ~ freedom, data = mixedEffect1)

#Replicating Figure 3
#NB: Figure 3 in the paper was generated from a slightly smaller number of elections in an earlier stage of the analysis. Because of that, the figure that is generated from the code above is minimally different. The substantive conclusions remain unchanged.

interactionvars <- c("pred","freedom01","education01_wi_01")
freedomdata <- data[,interactionvars]
freedomdata <- na.omit(freedomdata)

total <- nrow(freedomdata)
fr01 <- rep(0,total)

for (i in 1:total) {
if (freedomdata$freedom01[i] > .5) {fr01[i]=1}
}

fr01 <- factor(fr01, labels = c("Low Media Freedom", "High Media Freedom"))
freedomdata$fr01 <- fr01

edu <- ggplot(freedomdata, aes(x=freedomdata$education01_wi_01, y=freedomdata$pred))
edu <- edu + stat_smooth(method="lm", aes(linetype=fr01), shape=fr01, colour="black") + theme_bw() + scale_x_continuous("Education") + scale_y_continuous("Fitted Knowledge",breaks=c(.20,.21,.22,.23,.24,.25,.26,.27,.28,.29,.30,.31)) + theme(legend.title=element_blank(), legend.justification=c(0,.5), legend.position=c(0,.85), legend.background = element_rect())
edu <- edu + theme(axis.title.x = element_text(face="plain", color="black", size=15, hjust = .5,vjust = 0.1))
edu <- edu + theme(axis.title.y = element_text(face="plain", color="black", size=15, hjust = .5,vjust = 0.1))
edu <- edu + theme(legend.text = element_text(size=15))
print(edu)

#end script
